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Abstract -We investigate the optimization of synchronizability in multiplex networks and demon¬ 
strate that the interlayer coupling strength is the deciding factor for the efficiency of optimization. 
The optimized networks have homogeneity in the degree as well as in the betweenness central¬ 
ity. Additionally, the interlayer coupling strength crucially affects various properties of individual 
layers in the optimized multiplex networks. We provide an understanding to how the emerged 
network properties are shaped or affected when the evolution renders them better synchronizable. 
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. ^ Introduction: Synchronization is one of the most fas- 
^_cinating phenomena witnessed in interacting dynamical 
(-H units [I] . Given a wide applicability of this phenomenon, 
O i n the last two decades, a large number of studies have 
'“^teen devoted to understand the interplay of the individ- 
ual dynamical units and the underlying network struc- 
ture [5] , most importantly how structural properties affect 
the synchronizability dll]. Furthermore, various spectral 
properties of the underlying networks impart an under- 
j standing to the collective dynamical behavior [SHE], with 
the extremum eigenvalues of the corresponding Laplacian 
• matrices providing insight in to the synchronizability of a 
^network 

So far most of the works pertaining to the synchronization 
• • as well as optimally synchronized networks are restricted 
. ^ to the single layer networks, where attempts have been 
made to understand the evolutionary origin of the most 
synchronizable networks based on the minimization of the 
^ ratio (R) of the largest and the first nonzero eigenvalues 
of the Laplacian matrices using simulated annealing m- 
However, multilayer or multiplex networks have been in¬ 
creasingly realized to present a more realistic framework to 
model the interactions in real world systems m- Given 
the importance of the multilayer framework, there is a 
spurt in the activities of understanding and characteriz¬ 
ing their structural m and spectral properties [ISHUj. 
These works have shown that the synchronizability of mul¬ 
tilayer networks is dependent on the interlayer coupling 
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strength and attains a maximum value at some critical 
value of the strength, however, the degree of efficiency 
of optimization is not understood with respect to the in¬ 
terlayer coupling strength m- Furthermore, a weighted 
nature of the coupling strength has been realized in many 
real world networks mM- The optimization of synchro¬ 
nizability in weighted coupling is a complex problem due 
to the contribution arising from the degree of freedom in 
the coupling strength as well as the connection pattern of 
the networks. Few works pertaining to understand the im¬ 
pact of weighted nature of coupling strength on the syn¬ 
chronizability of a network conclude that the scale-free 
networks are more sensitive to various optimization meth¬ 
ods [inHin- 

In this Letter, we investigate the impact of structural 
properties of a network on its synchronizability under the 
multiplex network framework. The multiplex networks are 
evolved using the Genetic algorithm (GA) considering R 
as the fitness function. We find that the best synchroniz¬ 
able multilayer network has a stronger interlayer connec¬ 
tivity as compared to the connections within each layer. 
Further, the optimal network supports homogeneity in the 
degree, in addition to an emergence of the degree-degree 
correlations for the mirror nodes in different layers, re¬ 
flecting several distinguished features incorporating multi- 
plexity. We provide an understanding of various emerging 
structural features of the final evolved networks by relat¬ 
ing them with the load distribution on the interlayer links 
[5] as well as by the spectral properties of the underlying 
Laplacian matrices [MUSI- 
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Fig. 1: (a) Schematic diagram of a multiplex network con¬ 
sisting of two layers in which the solid lines represent in- 
tralayer couplings, whereas the dashed lines depict inter¬ 
layer couplings, (b) Schematic diagram denoting evolution 
of networks using the genetic algorithm, where Gi repre¬ 
sents the network in a population. G( and represent 
fitter parent and child networks, respectively. The solid 
arrow represents progression of generation from the initial 
population and the dashed arrow indicates that the fitter 
and child graphs form the initial population for the next 
generation. 


Theoretical framework: We consider a bi-layer mul¬ 
tiplex network, each layer having P network population 
with a fixed size N and a fixed average degree (fc). The 
networks representing different layers may have different 
architecture. Let the adjacency matrix of the first and sec¬ 
ond layers be denoted by A* and B®, respectively, where 
I ^ i ^ P. The weighted adjacency matrix M® of size 
2Nx2N for the multiplex network (Fig.[TJa)) is defined by 


M® 


A® D^I 
DJ B® 


( 1 ) 


where I is the identity matrix of size N x N and is the 
interlayer coupling strength. 

The Laplacian matrix T® obtained from IVP is defined by 


L®. = | ( 2 ) 

I —Mjf. otherwise, 

where d® = Let Aj < A 2 < .... ^ A^^ be the 

eigenvalues of L®. The eigenvalue ratio i?® = A^^/A^ of 
the T® matrix is used as a fitness function of multi¬ 
plex network in the genetic algorithm. A network having 
the lesser R value is more fit. Out of P population of 
the initial networks, we select P/2 fitter networks for the 
next generation, which we term as the parent networks 
(Fig. mb)). The remaining P/2 networks of the next gen¬ 
eration, known as the child networks, are constructed from 
the parent networks as follows. First a pair of the matri¬ 
ces is selected randomly from the P/2 population and then 
the adjacency matrices of these selected parent networks, 
termed as the parent matrices, are used to construct one 



Fig. 2: Minimized (circles) and initial (stars) average value 
of R for various interlayer coupling strength in panel (a) 
and its zoomed version in panel (b). Panel (c) shows the 
values of the correlation Lcorr for initial ER random net¬ 
works (stars) and optimized networks (circles). For each 
case, GA minimize the R till 6000 iterations and the av¬ 
erage is taken over population of 1000 networks and for 
the last 1000 iterations. Each initial layer of ER networks 
have (fc)=6 and N = 100 |26) . 


crossed child matrix. Each block of this child matrix is 
constructed from the block having the same position in 
one of the parent matrices which are selected with the 
equal probability. Then, by taking the upper triangular 
part of the crossed matrix, a symmetric child matrix is 
built leading to the undirected child network. Eurther, 
the average degree of the crossed matrix is preserved by 
randomly inserting or deleting 1 and 0 entries decided by 
a probability. This probability is dependent on the fluc¬ 
tuations of 1 entries from the adjacency matrices of the 
parent networks. The details of these steps are given in 
P5] . The process is repeated until we generate P/2 child 
networks which along with the P/2 parent networks form 
the same size of the population as of the initial networks. 
The networks are evolved until we attain the steady state. 
We measure various structural properties of networks as 
they evolve and analyze the interplay of these properties 
with the optimization of synchronizability. 


Results and Discussion: In order to draw a close coher¬ 
ence with the real world networks, we consider the sparse 
networks as the initial population in GA. Starting with 
networks having random architecture, we evolve them us¬ 
ing GA taking R as the fitness function. We find that at 
the lower interlayer couplings {Dx « 1), the R is high 
which decreases with an increase in the values. At 
a certain value, the R attains its minima, then sat¬ 
urates followed by an increase in its value for >> 1 
(Eig. [2Ka),(b)). Note that for all the simulations, the size 
and the average degree of individual layer remains same as 
for the initial network, hence, any change in the value of 
R for initial network population is arising due to changes 
in the values. The lower values indicate the higher 
values of R as at lower interlayer coupling strength the A 2 
values are low and at higher values the Xn values are 
high. Both the factors contribute to an increase in the R 
values m and for the model considered in the Letter, R 
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Fig. 3: Average value of Rnorm when individual network 
in a population is represented as (a) ER and (b) scale-free 
networks, respectively. For each case, GA minimizes the 
value of R till 6000 iterations and the average is taken over 
the population of 1000 networks for last 1000 iterations. 
Each layer of the networks have (fc)=6 and N = 100, which 
are preserved throughout the evolution. 



Fig. 4: Minimized (circles) and initial (stars) average value 
of R for various interlayer coupling strength in panel (a) 
and its zoomed version in panel (b). Panel (c) shows 
the values of the correlation Lcorr for initial SF networks 
(stars) and optimized networks (circles). For each case 
GA minimize the R till 6000 iterations and the average of 
is taken over population of 1000 networks and last 1000 
iterations. Each initial layer of scale-free networks have 
(fc)=6 and N = 100. 


can be determined as following: Case(l), when « 1 


R 


max \\Tnax{^R ) 4” Rx\ 
2Dx 


(3) 


where is the Laplacian of layer, \max{L°‘) is 
the maximum eigenvalue of the Laplacian of ath layer. 
Gase(2), when >> 1 


values of D^. The value of Rnorm starts increasing after 
attaining the minimum value. The above results indicate 
that the enhancement in the first term of the numerator 
of Eq. 131 which is proportional to the interlayer coupling 
strength, might lead to a decrease in the efficiency of the 
optimization (Rnorm)- The GA, thus maximizes the syn¬ 
chronizability for various values of the interlayer coupling 
strength (Fig. [SJa), (b)) with = 2 leading to the most 
efficient structure. 


^ 2D, + V2Xmax{L^^) ,,, 

- (■>) 

where is average Laplacian of two layers [16] . Eq. [3] 
indicates that i? is a decreasing function of for weaker 
interlayer coupling strength whereas Eq. |3] indicates that 
R is an increasing function of Dx for the stronger inter¬ 
layer coupling strength. Consequently there is a trade off 
between both the behaviors leading to a value of for 
which R achieves its minimum value. It turns out that 
this minimum value occurs at Dx = 2 and hence Eq. |3] 
stands as a better approximation for R. What follows 
that the structural changes in the architecture of the op¬ 
timized networks affect the second term in the numerator 
and the denominator which collectively contribute to min¬ 
imization of the R values. 

Further, we define the efficiency of optimization as 
Rnorm — Ropt/Rini^ where Rini and Ropt are the eigen¬ 
value ratio of the initial and the optimized networks, re¬ 
spectively. A lower value of Rnorm reflects a high ef¬ 
ficiency of the synchronizability. For the ER random 
networks, initially with an increase in the Dx values, 
Rnorm remains almost constant until the interlayer cou¬ 
pling strength reaches to 0.50 (Fig. [3[a)). However, the 
value of Rnorm being less than one implies that there is 
an optimization efficiency with respect to the initial pop¬ 
ulation. Thereafter, it exhibits an abrupt decrease and at¬ 
tains a minimum value for Dx = 2. Therefore, GA finds a 
variability in the efficiency of the optimization for various 


In order to capture the structural features emerging in the 
best synchronizable networks, we define the Pearson prod¬ 
uct correlation (Lcorr) between the degrees of the mirror 
nodes in the pair of the layers as the evolution progresses 
as; 

r (S.(fc.-(fc))(S.(^.-(A0) 

where ki and Ki are the degrees of node in the first 
and the second layers, respectively. The terms (fe) and 
(K) denote the average degrees of the first and the sec¬ 
ond layers, respectively. Lcorr is the average over the Lcorr 
values of the population used in the GA. For lower val¬ 
ues of Dx{< 0.5), the optimized networks do not reveal 
significant changes in Lcorr and remain close to zero, with 
small random fluctuations. With a further increase in Dx 
values, the average correlation Lcorr decreases for the op¬ 
timized networks (Fig. [2Ic)). This can be considered to 
occur in order to minimize the value of maximal load over 
the interlayer links. 

The reasons behind the emerging behavior of Lcorr might 
be explained as follows. When the higher degree nodes in 
one layer are connected with the higher degree nodes of the 
second layer and the interlayer coupling strength is high 
as well {Dx >1), there is a very high load on the interlayer 
links, making it to become a dominant factor for affecting 
the synchronizability of the network |5]. Whereas, in a sit¬ 
uation when the higher degrees in one layer are connected 
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Fig. 5: (a) Average values of Rnorm as a function of D^. 
For each case GA minimize the Rnorm till 6000 iterations 
and the average of is taken over population of 1000 net¬ 
works and last 1000 iterations. Panel (b) shows the con¬ 
vergence of average of Rnorm values over GA networks 
populations with generations(G) for = 1. For both 
cases Ey = A and each initial layer of ER networks have 
(k)=6 and N = 100. 


with the lower degree nodes in another layer, the value 
of the maximal load on the interlayer links becomes less, 
leading to a homogeneity of the loads on the interlayer 
links. In order to achieve this homogeneity, at higher val¬ 
ues of Dx , the optimized networks exhibit negative values 
of Acorr ■ These results demonstrate an impact of the mul- 
tiplexity on the maximization of the synchronizability as 
well as an impact of the weighted couplings on the evolved 
network architecture. 

Further, in order to demonstrate the impact of structural 
properties of the initially selected networks population on 
the structural properties of the optimized networks, we 
present the results for the scale-free networks taken as 
the initial population. The scale-free networks are gener¬ 
ated using Barabasi-Albert algorithm [53]. As observed 
for the ER networks, the GA approach successfully mini¬ 
mizes the R values for various values for the scale-free 
networks as well. The value of D^, for which R takes a 
minimum value, also remains same (Fig. jT^a), (b)). The 
Lcorr value for the initial networks population is positive 
(Fig. HKc)) since the mirror nodes in both the layers are 
associated with the same time in the growth algorithm, 
but the optimized networks exhibit negative L^orr values 
as Dx increases (Fig. |3Kc)). The efficiency of optimization 
for scale-free networks taken as the initial population, cap¬ 
tured through Rnorm, also exhibits a behavior similar to 
the ER random networks (Fig.|2Kb)) at the same interlayer 
coupling strength, though the values of Rnorm are lesser 
than those of ER random networks probably due to the 
higher values of R for scale-free networks as compared to 
the ER random networks. 

In order to analyze the impact of intralayer coupling 
strength on the global synchronizability, we introduce a 
parameter Ey. The new weighted adjacency matrix A'R 



Fig. 6: Average value of the variance of degree and be¬ 
tweenness centrality of the nodes in the optimized net¬ 
works in the left and right panels, respectively. Initial 
layer of each ER network has {k)=6 and N = 100. 


for the multiplex network can be defined as. 




A* Dxl 
Dxl EyB'- 


(5) 


We consider one layer having intralayer coupling strength 
stronger than that of the other layer {Ey = 4 in Eq. |5|). 
With an initial increase in D^, the values Rnorm remains 
unaffected. After a certain value of D^, Rnorm decreases 
drastically and thereafter for a certain range of the 
rate of decrement becomes less. After attaining a mini¬ 
mum value, Rnorm increases again (Fig. 5(a)). This be¬ 
havior shares a close similarity to the one observed for Ey 
= 1 (Eq. 3), the major difference being the shift in the 
minima of Rnorm towards right, appearing at = 16. 
The probable reason behind this shift might be rescaling 
of the Dx values with respect to Ey towards the regime 
Dx « 1. As a result, Eq. 3 holds good instead of Eq. 4 
for the rescaled values of Dx lying in the regime Dx << 
1 . 

Furthermore, there is a nearly similar behavior of conver¬ 
gence of Rnorm with an increase in generation (G) of the 
GA. We find that initially with an increase in the gener¬ 
ation, the Rnorm values exhibit a faster decrement (Fig. 
5(b)). After a certain value of G, the Rnorm converges 
to a hxed value and the properties affecting the optimiza¬ 
tion of the networks population saturate. This leads to 
a steady value of Rnorm which remains unaffected by a 
further increase in G. 

Further, to analyze the properties in the optimized net¬ 
works, we calculate the variance in the degree and the 
betweenness centrality [5S] for the different Dx values. 
The initial ER random networks population has cer¬ 
tain amount of degree heterogeneity due to the Poisson 
degree distribution [24] . For the lower interlayer cou¬ 
pling strength {Dx < 2), the heterogeneity is suppressed 
(Fig. |n|). With a further decrease in the Dx values, the 
suppression in the degree heterogeneity remains invariant, 
beyond which the suppression decreases in the optimized 
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networks (Fig. EKa)). Note that there is a higher hetero¬ 
geneity in the connectivity for = 2 as compared to 
= 1. Further, the betweenness centrality of a node can also 
be treated as a measure of the load distribution over the 
nodes in a network |12j . The ER random networks forming 
the initial population again display a certain amount of the 
betweenness centrality heterogeneity. For the higher cou¬ 
pling strength, the heterogeneity in betweenness centrality 
is not much suppressed in the optimized networks, how¬ 
ever with a decrease in the interlayer coupling strength the 
suppression increases and converges to the homogeneous 
distribution at the lower coupling regime (Fig. EKb)). 
Though the optimized networks suppress heterogeneity in 
the degree as well as in the betweenness centrality, the 
interlayer coupling strength plays an instrumental role in 
deciding the amount of heterogeneity. The larger inter¬ 
layer coupling strengths lead to more heterogeneity in the 
various structural properties of the optimized networks. 

These behaviors can be explained from the nature of Eq. [3] 
and Eq. |4l From Eq. [3j which holds good in the regime 
Dx >> 1, a structural change leads to a variation in the 
maxima of \max over layers by keeping its denominator 
constant for a given Dx value. Since the upper bound of 
^max is an increasing function of the maximum degree of a 
network |27] , a minimized value of Xmax is achieved when 
the maxima of the degree of networks in both the layers 
are minimized. However in regime Dx >> 1, average value 
of A 2 of the Laplacian for the layers appear in the denomi¬ 
nator of Eq.m which is lower bounded by the inverse of the 
diameter of the network multiplied with the total number 
of connections [19], providing another degree of freedom. 
As a result sufficient amount of the homogeneity does not 
emerge in the optimized networks. 

Conclusion: To conclude, we study the behavior of syn¬ 
chronizability in the multiplex networks with respect to 
interlayer coupling strength. Various bounds for Dx pro¬ 
vide an understanding to how interlayer coupling strength 
turns out to be a deciding factor for achieving the most 
synchronizable networks as well as for the efficiency of the 
optimization. The results presented here have two fold 
applications, one being the influence of the weights in cou¬ 
pling strength on the synchronizability and the other per¬ 
tains to the importance of the multiplex network frame¬ 
work in capturing the enhancement in the synchronizabil¬ 
ity and the identification of the structural features whose 
interplay can help in designing the best synchronizable 
networks. The behavior of the synchronizability with re¬ 
spect to the interlayer coupling strength is robust against 
changes in the initial networks architecture. Based on the 
eigenvalue ratio, we argue that a particular interlayer cou¬ 
pling strength value leads to the highest efficiency of the 
optimization. Interestingly, the dynamical behavior of in¬ 
dividual layers has been demonstrated to be largely af¬ 
fected while using them in the multiplex network frame¬ 
work. While the individual layers may not be favoring 
synchronizability, the evolved multiplex networks can be 


best synchronizable demonstrating the importance of mul¬ 
tilayer framework. Additionally, the betweenness central¬ 
ity of the edges being treated as the load distribution m 
explains the emergence of the negative correlation in the 
mirror nodes in the best synchronizable networks. 

The framework and analysis presented here can be ex¬ 
tended further in identifying a particular layer of the mul¬ 
tiplex network which can be targeted to achieve the best 
or the worst synchronizable multiplex networks. Eurther, 
the approach presented in this Letter can be adopted to 
optimize other dynamical properties [28] in multiplex net¬ 
works having restrictions on wiring topologies |29] . 
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